//    Copyright (C) Kherty.  All rights reserved.
using System;

namespace OpenLS.Spreadsheet.StandardAddIn.Statistics
{
    internal delegate double GetFn(int n, double x);

    internal class ContinuedFraction
    {
        private static readonly double DefaultEpsilon = 10e-9;

        private readonly GetFn getA;
        private readonly GetFn getB;

        public ContinuedFraction(GetFn getA, GetFn getB)
        {
            this.getA = getA;
            this.getB = getB;
        }

        public double Evaluate(double x)
        {
            return Evaluate(x, DefaultEpsilon, int.MaxValue);
        }

        public double Evaluate(double x, double epsilon)
        {
            return Evaluate(x, epsilon, int.MaxValue);

        }

        public double evaluate(double x, int maxIterations)
        {
            return Evaluate(x, DefaultEpsilon, maxIterations);
        }

        public double Evaluate(double x, double epsilon, int maxIterations)
        {
            double p0 = 1.0;
            double p1 = getA(0, x);
            double q0 = 0.0;
            double q1 = 1.0;
            double c = p1/q1;
            int n = 0;
            double relativeError = double.MaxValue;
            while (n < maxIterations && relativeError > epsilon)
            {
                ++n;
                double a = getA(n, x);
                double b = getB(n, x);
                double p2 = a*p1 + b*p0;
                double q2 = a*q1 + b*q0;
                if (double.IsInfinity(p2) || double.IsInfinity(q2))
                {
                    // need to scale
                    if (a != 0.0)
                    {
                        p2 = p1 + (b/a*p0);
                        q2 = q1 + (b/a*q0);
                    }
                    else if (b != 0)
                    {
                        p2 = (a/b*p1) + p0;
                        q2 = (a/b*q1) + q0;
                    }
                    else
                    {
                        // can not scale an convergent is unbounded.
                        throw new ConvergenceException();
                    }
                }
                double r = p2/q2;
                relativeError = Math.Abs(r/c - 1.0);

                // prepare for next iteration
                c = p2/q2;
                p0 = p1;
                p1 = p2;
                q0 = q1;
                q1 = q2;
            }

            if (n >= maxIterations)
            {
                throw new ConvergenceException();
            }
            return c;
        }
    }
}